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LONG-TERM  GOALS 

The  long-term  goal  of  this  research  is  to  construct  global  and  mesoscale  nonhydrostatic  numerical 
weather  prediction  (NWP)  models  for  the  U.S.  Navy  using  new  numerical  methods  specifically 
designed  for  modem  computer  architectures.  To  take  full  advantage  of  distributed-memory  computers, 
the  global  domains  of  these  new  models  are  partitioned  into  local  sub-domains,  or  elements,  which 
can  then  be  solved  independently  on  multiple  processors.  The  numerical  methods  used  on  these  sub- 
domains  are  local,  high-order  accurate,  fully  conservative,  and  highly  efficient.  Using  these  ideas  we 
are  developing  global  and  mesoscale  nonhydrostatic  atmospheric  models  that  will  improve  upon  the 
operational  models  currently  used  by  all  U.S.  agencies  including  the  U.S.  Navy. 

OBJECTIVES 

The  objective  of  this  project  is  to  construct  new  high-order  local  methods  for  the  Navy’s  next- 
generation  global  and  mesoscale  nonhydrostatic  NWP  models.  The  high-order  accuracy  of  these 
methods  will  ensure  that  the  new  model  yields  better  forecasts  than  the  current  global  spherical 
harmonics  model  (NOGAPS)  and  better  accuracy  than  the  current  mesoscale  finite  difference  model 
(COAMPS).  The  objective  is  to  achieve  this  accuracy  while  increasing  the  geometric  flexibility  to  use 
any  grid  as  well  as  to  increase  the  efficiency  of  these  models  on  large  processor-count  distributed- 
memory  computers.  Higher  efficiency  means  that  the  new  models  will  require  less  computing  time 
which  then  allows  for  increasing  the  number  of  ensemble  members  and/or  increasing  the  resolutions 
of  the  NWP  models.  The  methods  that  we  propose  to  use  for  these  models  are  state-of-the-art  and  are 
not  being  used  by  either  current  or  newly  emerging  NWP  and  climate  models. 

APPROACH 

To  meet  our  objectives  we  explore: 

1.  spectral  element  (SE)  and  discontinuous  Galerkin  (DG)  spatial  discretization  methods; 

2.  high-order  semi-implicit  (SI)  time-integrators  with  adaptive  time-stepping  for  vastly  improved 
efficiency; 

3.  high-order  Lagrangian-like  time-integrators  that  are  fully  conserving  and  scale  well  on  modem 
computer  architectures; 
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Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 
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4.  various  forms  of  the  governing  equations  in  order  to  maximize  accuracy,  efficiency,  stability, 
and  conservation  properties; 

5.  unified  hydrostatic/nonhydrostatic  formulations  of  the  equations;  and 

6.  fully  unstructured  (and  possibly  adaptive)  grids. 

The  power  of  SE  and  DG  methods  is  that  they  are  high-order  accurate  yet  are  completely  local  in 
nature  -  meaning  that  the  equations  are  solved  independently  within  each  individual  element  and 
processor.  Furthermore,  high-order  methods  have  minimal  dispersion  error  -  this  is  an  important 
property  for  capturing  fine-scale  atmospheric  phenomena  (e.g.,  tropical  cyclones,  Kelvin  and  Rossby 
waves).  The  theoretical  development  of  SE  and  DG  methods  are  now  well-established  and  these 
methods  are  currently  the  two  most  successful  methods  found  in  the  literature  for  fluid  flow  problems. 

Semi-implicit  (SI)  and  Lagrangian  time-integrators  offer  vast  improvements  in  efficiency  due  to  the 
longer  time  steps  that  they  permit;  it  should  be  mentioned  that  semi-implicit  and  Lagrangian-like 
methods  can  be  classified  together  under  the  heading  of  implicit-explicit  (IMEX)  methods  which  has 
garnered  much  attention  in  the  computational  mathematics  literature.  Furthermore,  in  order  to  reap  the 
full  benefits  of  the  high-order  spatial  discretization  methods  requires  increasing  the  order  of  accuracy 
of  the  time-integration  methods  as  well;  this  is  a  topic  that  too  often  has  been  ignored  by  most 
scientific  computing  communities,  including  the  NWP  community.  Lagrangian  methods  have  not  been 
used  successfully  for  mesoscale  modeling  because  of  their  lack  of  conservation.  Another  problem  that 
they  pose  is  that  they  require  vast  amounts  of  inter-processor  communication  on  a  distributed-memory 
computer.  We  have  worked  on  Lagrangian-like  methods  that  are  conserving  and  require  no  additional 
inter-processor  communication  (see  [11]). 

Before  committing  resources  towards  the  development  of  new  NWP  models,  it  is  important  to  identify 
the  form  of  the  governing  equations  that  is  most  capable  of  conserving  all  quantities  deemed 
important.  We  performed  a  study  on  this  topic  this  year-  that  is,  to  identify  the  form  of  the  governing 
equations  capable  of  representing  conservation  of  either  mass,  energy,  or  both.  In  addition,  we  have 
begun  to  analyze  various  forms  of  the  governing  equations  with  respect  to  robustness,  flexibility,  and 
efficiency  in  the  context  of  implicit-explicit  (IMEX)  time-integration  methods.  Within  this  work  we 
have  also  explored  hybrid  models  that  solve  either  the  hydrostatic  or  nonhydrostatic  equations.  This 
feature  allows  the  models  to  be  used  for  research  purposes  by  Navy  scientists  in  order  to  test  the 
importance  of  nonhydrostatic  phenomena  at  specific  resolutions. 

One  final  area  that  needs  to  be  explored  is  the  concept  of  adaptive  grids.  In  the  past  few  years,  adaptive 
grids  have  gained  considerable  momentum  in  the  atmospheric  modeling  community  -  in  fact,  I  have 
been  invited  to  give  a  keynote  lecture  at  the  University  of  Reading  in  March  2009  to  kick-off  a  year 
long  program  on  adaptive  modeling  at  the  Newton  Institute  in  Cambridge  University,  England. 

WORK  COMPLETED 

In  this  section,  we  describe  the  work  completed  this  fiscal  year.  The  work  can  be  categorized  into  three 
sections:  global  modeling,  mesoscale  modeling,  and  development  of  new  numerical  machinery  for 
solving  mathematical  issues  related  to  both  global  and  mesoscale  modeling. 
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Global  Modeling  Work.  Work  on  the  global  model  consisted  of  incremental  improvements  to  the 
model  that  has  been  under  development  for  the  past  few  years.  This  model,  NSEAM,  is  hydrostatic  and 
therefore  would  have  a  relatively  short  life-span  as  an  operational  model  since  we  are  developing 
nonhydrostatic  models  for  both  global  and  mesoscale  applications.  Nonetheless,  NSEAM  has  been 
shown  to  be  an  effective  tool  for  studying  the  strengths  and  weaknesses  of  the  NOGAPS  physical 
parameterization  package.  Using  NSEAM,  we  have  published  a  paper  (see  ref.  [6])  showing  the 
sensitivity  of  Kelvin  waves  and  the  Madden- Julian  Oscillation  on  the  NOGAPS  physics;  this  paper 
exposed  the  strong  sensitivity  of  the  current  physical  parameterizations  on  the  vertical  discretization; 
this  includes  sensitivity  to  the  resolution  as  well  as  to  the  placement  of  the  grid  points.  This  strong 
sensitivity  to  the  vertical  structure  of  global  models  is  an  issue  that  needs  to  be  addressed  by  the 
atmospheric  modeling  community. 

The  path  to  building  a  nonhydrostatic  global  model  can  take  one  of  two  paths:  the  first  is  to  extend  the 
existing  hydrostatic  NSEAM  model  to  be  nonhydrostatic  and  the  second  approach  is  to  build  a  state- 
of-the-art  nonhydostatic  mesoscale  model  that  can  then  be  modified  to  run  as  a  global  model.  The  first 
approach  would  constrain  the  equations  to  be  written  in  a  form  that  would  not  be  able  to  conserve  mass 
formally.  This  is  the  reason  we  have  decided  on  the  second  approach.  Let  us  now  describe  the  work 
completed  this  year  on  the  mesoscale  models. 

Mesoscale  Modeling  Work.  The  We  have  continued  our  analysis  of  the  various  forms  of  the  Euler 
equations  and  their  advantages/disadvantages  for  mesoscale  modeling.  We  have  included  a  few  more 
forms  to  our  study  totaling  5  different  forms  of  the  equations  which  we  now  discuss.  Specifically  we 
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analyzed  the  following  equations:  Set  1  is  defined  as  follows  where  the  solution  vector  is  Exner 
pressure,  velocity,  and  potential  temperature.  Set  1  is  the  equation  set  used  in  the  U.S.  Navy’s 
mesoscale  model  COAMPS.  The  main  problem  with  this  equation  set  is  that  it  cannot  conserver  either 
mass  or  energy. 


Set  2  in  conservation  form  (denoted  as  Set  2C)  is  defined  as  follows  where  the  solution  vector  is 
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density,  momentum,  and  density  potential  temperature.  Set  2  is  the  form  used  in  WRF.  This  form  is 
very  attractive  because  it  conserves  mass  although  it  does  not  conserve  energy.  It  does,  however, 
conserve  density  potential  temperature  which  is  related  to  entropy. 


This  set  is  of  interest  because  it  can  also  be  written  in  non-conservation  form  while  still  conserving 
mass.  We  shall  refer  to  it  as  Set  2NC  (for  non-conserving),  and  it  is  defined  as  follows: 
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The  interest  in  equation  Set  2NC  is  that  it  conserves  mass  and  offers  much  flexibility  in  the  type  of 
time-integrators  that  can  be  used  with  it.  For  example,  note  that  the  first  two  terms  (in  red  font)  of  each 
of  the  components  of  mass,  momentum,  and  potential  temperature  can  be  recast  as  a  Lagrangian 
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derivative.  This  then  allows  the  use  of  Lagrangian-like  time-integrators  to  be  used  such  as  those  we 
describe  in  [11]. 


Set  3  is  defined  as 
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where  the  solution  vector  is  density,  momentum,  and  density  total  energy,  where  E=pe  with  e=cvT 
+0.5u*u  +  (|),  in  other  words  e  represents  internal  energy,  kinetic  energy,  and  potential  energy.  This 
equation  set  is  not  used  in  atmospheric  modeling  but  is  the  equation  of  choice  in  computational  fluid 
dynamics  (CFD).  This  set  is  very  attractive  because  it  conserves  both  mass  and  energy  regardless  of 
whether  the  flow  is  inviscid  or  viscous  (with  the  proper  viscous  stressed  included)  (see  [9]).  One 
question  we  had  about  this  set,  however,  is  whether  it  could  be  coupled  to  existing  physical 
parameterization  packages  which  rely  on  potential  temperature  and  not  energy;  this  question  has  been 
answered  in  the  past  year  with  the  recent  results  of  the  Japanese  N1CAM  model  which  in  fact  uses  this 
equation  set. 

The  final  equation  set  studied  is  Set  4  which  is  written  as  follows: 

—  +  u*  V/7+  pV  •  u  =  0;  —  +  u*Vu  +  —  VP  =  -/(kxu)-gk;  —  +  u  •  VP  +  yPV  •  u  =  0 
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Note  that  this  equation  is  also  amenable  to  various  time-integration  methods  including  Lagrangian 

methods  since  the  first  two  terms  of  each  of  the  components  of  the  equations  can  be  written  as  a 

Lagrangian  derivative.  Furthermore,  as  in  Sets  2C  and  2NC,  Set  4  can  formally  conserve  mass  but  not 

energy.  Set  4  is  a  very  good  compromise  between  conservation  and  efficiency.  Note  that  Sets  1  and  4 

are  the  only  fully  closed  systems  requiring  no  equation  of  state,  whereas  Sets  2C,  2NC,  and  3  all 

require  an  equation  of  state  in  order  to  couple  the  extra  variable  (pressure)  to  the  prognostic  variables. 

Set  4  is  the  equation  set  used  by  the  LM  model  of  the  German  Weather  Service  and  has  some  good 

properties  that  we  are  analyzing. 

Using  these  5  equation  sets  we  developed  x-z  slice  mesoscale  models  to:  compare  the  spectral  element 
and  discontinuous  Galerkin  methods,  analyze  semi-implicit  time-integrators,  and  to  see  how  these 
models  behaved  under  a  series  of  test  cases  including  sharp  front  simulations  and  nonhydrostatic  flow 
over  mountains.  This  work  resulted  in  a  peer-reviewed  article  that  appeared  this  year  (see  ref.  [9])  and 
another  that  is  in  preparation. 

Mathematical  Issues  for  NWP  Modeling.  In  collaboration  with  Marco  Restelli  (a  previously  funded 
ONR  Global  VSP  visitor  at  NPS)  the  PI  (Giraldo)  was  able  to  develop  the  first  semi-implicit  time- 
integrator  for  any  DG  model.  This  result  is  very  important  because,  until  this  point,  DG  models  were 
not  competitive  with  other  gridpoint  models  in  terms  of  efficiency  because  it  was  not  known  how  to 
construct  semi-implicit  time-integrators  for  DG.  This  work  is  currently  under  review  (see  ref.  [5]). 

The  high-order  in  time  generalization  of  this  method  was  submitted  to  a  journal  this  year  (see  ref.  [1]). 
In  work  in  progress,  we  are  extending  this  method  to  handle  hydrostatic/nonhydrostatic  hybrid 
equations  (that  allow  the  study  to  determine  when  nonhydrostatic  effects  are  important)  and  adaptive 
time-stepping  (that  allow  the  model  to  automatically  modify  its  time-step  size  in  order  to  remain  stable 
under  any  flow  situation  and  yield  efficient  time-integrations).  This  work  is  in  collaboration  with 
another  ONR  Global  VSP  visitor,  Matthias  Lauter. 
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Another  interesting  mathematical  issue  for  NWP  modeling  that  we  explored  this  year  involved  the 
development  of  high-order  non-reflective  boundary  conditions  for  mesoscale  models.  This  work  was 
conducted  by  Air  Force  Maj.  John  Dea  for  his  PhD  dissertation  which  he  successfully  defended  in 
September  2008  under  my  supervision.  The  results  show  that  using  a  characteristic  decomposition  of 
the  governing  equations  allows  for  a  relatively  straightforward  implementation  of  Higdon-type  non- 
reflecting  boundary  conditions  (NRBCs)  for  the  Euler  equations;  Higdon  NRBCs  are  the  high-order 
generalization  of  the  classical  Sommerfeld  conditions.  This  work  has  been  submitted  to  peer-reviewed 
journals  (see  refs.  [2,3,8]).  In  his  dissertation,  Dr.  Dea  developed  the  theoretical  framework  for 
implementing  this  idea  and  implemented  it  to  the  linearized  Euler  equations  (for  stratified  flow 
relevant  to  the  atmosphere)  using  the  finite  difference  method.  A  new  PhD  student  will  now  extend  Dr. 
Dea’s  ideas  to  high-order  spatial  discretization  and  time-integration  methods  in  order  to  test  them  on 
our  new  x-z  nonhydrostatic  mesoscale  models.  In  addition,  a  National  Research  Council  (NRC)  post¬ 
doc  (Dr.  Jim  Kelly)  will  work  on  a  number  of  mathematical  topics  related  to  the  mesoscale  models 
including:  NRBCs,  optimal  time-integrators,  parallelization  of  the  models,  and  extending  them  to  three 
dimensions. 

RESULTS 

Mesoscale  Nonhydrostatic  Atmospheric  Models.  In  Ref.  [9]  we  present  five  mesoscale  models  using 
equation  sets  1,  2,  and  3,  in  conjunction  with  the  SE  and  DG  methods.  For  example,  SE1  means  that 
the  SE  method  was  used  for  set  1.  Note  that  there  is  no  DG1  because  set  1  is  not  in  conservation  form 
and,  therefore,  cannot  be  used  with  the  DG  method.  In  Ref.  [9]  we  show  results  for  six  test  cases 
ranging  from  rising  bubbles,  density  currents,  to  mountain  wave  problems.  In  Ref.  [5]  we  show  results 
for  the  DG3  with  the  semi-implicit  method.  We  now  show  results  from  both  of  these  papers.  Figure  1 
shows  the  results  for  a  cold  bubble  dropping  in  an  isothermal  atmosphere.  Figure  la  (left  panel)  shows 
the  potential  temperature  perturbation  contours  for  4  meter  resolution  after  100  seconds  and  Figure  lb 
(right  panel)  shows  the  results  after  200  seconds.  The  results  show  that  the  model  resolves  the  fine 
scale  features  extremely  well,  even  when  using  such  large  time-steps  -  moreover,  the  results  are  fully 
conservative  meaning  that  both  mass  and  energy  are  conserved  up  to  machine  precision. 
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a)  b) 

Figure  1:  Potential  Temperature  contours  for  a  cold  bubble  test  after  a)  100  seconds  and  b)  200 
seconds  for  the  semi-implicit  DG3  model  with  a  Courant  number  of  100. 


In  Fig.  2  we  show  the  performance  study  of  the  semi-implicit  discontinuous  Galerkin  method  for  the 
cold  bubble  test.  In  order  to  construct  an  efficient  semi-implicit  formulation,  after  the  system  of 
equations  are  discretized  in  time,  they  must  be  reduced  to  a  block  LU  decomposition  for  one  of  the 
variables  only.  The  resulting  equation  takes  the  form  of  a  Helmholtz  operator  and  for  this  reason,  this 
reduction  is  called  the  pseudo-Helmholtz  form  of  the  equations.  Figure  shows  the  number  of  GMRES 
iterations  per  time-step  (•  on  left-axis),  the  total  number  of  GMRES  iterations  (  ♦  on  first  right-axis), 
and  total  CPU  time  (*  on  second  right-axis)  as  functions  of  Courant  number.  Figure  2a  (left  panel) 
shows  the  results  using  the  pseudo-Helmholtz  form  while  Figure  2b  (right  panel)  shows  the  results 
without  the  pseudo-Helmholtz  form.  The  main  point  of  these  results  is  that  the  total  CPU  time  does 
not  grow  linearly  with  Courant  number  if  one  uses  the  pseudo-Helmholtz  form,  without  the  pseudo- 
Helmholtz  form,  the  cost  clearly  increases  linearly  with  increasing  Courant  number  (i.e.,  time-step 
size). 


a)  b) 

Figure  2:  Performance  study  of  the  semi-implicit  DG3  model  for  the  cold  bubble  test  for  various 
Courant  numbers  a)  with  and  b)  without  the  pseudo-Helmholtz  form. 
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Figure  3  shows  the  vertical  velocity  contours  and  normalized  momentum  flux  for  the  linear 
hydrostatic  mountain  test.  In  Fig.  3a  (left  panel)  the  dashed  red  lines  denote  the  analytic  solutions 
while  the  solid  black  lines  are  the  numerical  solutions.  Note  how  accurately  the  numerical  results 
agree  with  the  analytic  solutions.  In  Fig.  3b  (right  panel)  we  show  the  normalized  momentum  flux 
which  is  a  good  measure  of  the  accuracy  of  the  model.  A  perfect  model  would  yield  values  of  1  and 
we  can  see  that  we  are  near  1  for  all  five  models  developed  and  for  the  entire  vertical  column  of  the 
atmosphere.  The  values  near  12  kilometers  differ  from  the  value  1  due  to  the  effects  of  the  non- 
reflecting  boundary  conditions. 


a)  b) 

Figure  3:  Results  for  the  linear  hydrostatic  mountain  test.  We  show  a)  vertical  velocity  contours 
showing  the  analytic  solution  (dashed  red  lines)  and  the  numerical  solution  (solid  black  lines)  and 
b)  the  normalized  momentum  flux  for  5  different  models  we  have  developed. 


Figure  4  shows  the  vertical  velocity  contours  and  normalized  momentum  flux  for  the  linear 
nonhydrostatic  mountain  test.  In  Fig.  4a  (left  panel)  the  dashed  red  lines  denote  the  analytic  solutions 
while  the  solid  black  lines  are  the  numerical  solutions.  Note  how  accurately  the  numerical  results  once 
again  agree  with  the  analytic  solutions.  In  Fig.  4b  (right  panel)  we  show  the  normalized  momentum 
flux  for  all  five  models  developed  and  we  see  that  we  do  indeed  get  values  near  1 . 


a)  b) 

Figure  4:  Results  for  the  linear  nonhydrostatic  mountain  test.  We  show  a)  vertical  velocity  contours 
showing  the  analytic  solution  (dashed  red  lines)  and  the  numerical  solution  (solid  black  lines)  and 
b)  the  normalized  momentum  flux  for  5  different  models  we  have  developed. 
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From  the  preliminary  results  of  the  previous  year,  we  concluded  that  Sets  2  and  3  are  clearly  superior 
to  Set  1.  However,  this  year  we  have  included  two  additional  sets  (Set  2NC  and  Set  4).  The  reason  for 
considering  these  two  sets  is  to  consider  not  only  the  accuracy  of  the  equation  sets  with  respect  to  the 
SE  and  DG  spatial  discretization  but  also  to  consider  the  effects  of  the  time-integrators.  While  the 
preliminary  results  of  this  year  show  that  all  equation  sets  yield  similar  results  (in  terms  of  accuracy, 
efficiency,  and  robustness)  set  2NC  and  set  4  offer  the  most  flexibility  in  the  choices  of  time- 
integrators  that  can  be  used  with  them;  we  plan  to  demonstrate  this  in  the  next  year.  Furthermore,  we 
propose  to  develop  Eulerian  and  Lagrangian  semi-implicit  time-integrators  for  these  equation  sets. 

IMPACT/APPLICATIONS 

NOGAPS  and  COAMPS  are  run  operationally  by  FNMOC  and  is  the  heart  of  the  Navy’s  operational 
support  to  nearly  all  DOD  users  worldwide.  This  work  targets  the  next-generation  of  these  systems  for 
massively  parallel  computer  architectures.  NSEAM  and  its  mesoscale  cousins  have  been  designed 
specifically  for  these  types  of  computer  architectures  while  offering  more  flexibility,  robustness,  and 
accuracy  than  the  current  operational  systems.  Additionally,  the  new  models  are  expected  to  conserve 
all  quantities  such  as  mass  and  energy  and  use  state-of-the-art  time-integration  methods  that  will 
greatly  improve  the  capabilities  of  the  Navy’s  forecast  systems. 

TRANSITIONS 

Improved  algorithms  for  model  processes  will  be  transitioned  to  6.4  (PE  0603207N)  as  they  are  ready, 
and  will  ultimately  be  transitioned  to  FNMOC  with  future  NOGAPS  upgrades. 

RELATED  PROJECTS 

Some  of  the  technology  developed  for  this  project  could  be  used  to  improve  NOGAPS  in  other  NRL 
projects.  The  work  performed  in  this  work  unit  on  coupling  NSEAM  with  NOGAPS  physics  has 
already  revealed  some  sensitivities  of  the  physical  parameterization  to  the  vertical  coordinate;  this 
information  can  now  be  used  to  improve  the  forecasts  of  NOGAPS.  In  addition,  the  work  on  the 
mesoscale  models  will  help  improve  CO  AMPS.  An  example  is  the  time-integration  methods  that  we 
are  exploring  for  the  new  models  may  well  be  incorporated  into  the  current  operational  version  of 
COAMPS. 
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